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Abstract 

In the hadronic phase, the dominant configuration of QCD with two flavors of massless quarks 
is a gas of massless pions. We calculate the bulk viscosity (C) using the Boltzmann equation with 
the kinetic theory generalized to incorporate the trace anomaly. We find that the dimensionless 
ratio C/s, s being the entropy density, is monotonic increasing below T =120 MeV, where chiral 
perturbation theory is applicable. This, combined with previous results, shows that C,/ s reaches 
its maximum near the phase transition temperature T^., while r]/s, rj being the shear viscosity, 
reaches its minimum near Tc in QCD with massless quarks. 



I. INTRODUCTION 



Transport coefficients of Quantum Chromodynamics (QCD) are of high interests re- 
cently. This was triggered by the discovery that quark gluon plasma (QGP) has a viscosity 
close to the conjectured universal minimum bound jj, indicating that QGP is close to 
a "perfect fluid" , |j| just above the deconfinement temperature. This bound, rj/s 
> l/47r, s being the entropy density, is motivated by the uncertainty principle and is 
found to be saturated for a large class of strongly interacting quantum field theories whose 
dual descriptions in string theory involve black holes in anti-de Sitter space 
lere are some debates about whether the minimum bound on rj/s is truly universal 



ll| and the RHIC data might be better fit with rj/s < l/Air |l2l . (lattice results 



for gluon plasma, however, is still consistent with the bound [Ij, Il5|). In any case, smaller 
rj implies stronger interparticle interaction (here r] is normalized by the density) and the 
smallness of QGP t] indicating an intriguing strongly interacting state is reached near the 
deconfinement temperature. 

In general, the minimum of r?/s is found near the phase transition temuerature Tn or 

nnn 

when the system goes through a fast crossover. This behavior was observed Ha, 117[ in all 
the materials, including N, He, and II2O, with data available in the NIST and CODATA 



websites ilA uM- Surprisingly, it is also observed in QCD at zero chemical potential 

nn nn 

near the nuclear liquid gas phase transition [16|, |21[ , and in cold fermionic atom 



systems at the limit with two-body scattering length tuned to infinity [22j]. Using weakly- 
coupled real scalar field theories, in which perturbation is reliable, the same rj/s behaviors 
in first-, second-order phase transitions and crossover also emerge as in the liquid-gas 
transitions in N, He, and H2O and essentially all the matters with data available in the 



NIST database mentioned above 



23|. This agreement is expected to hold when the theory 



is generalized to components with an 0{N) symmetry. Thus, these behaviors might 
be general properties of fluid and might be used to probe the QCD critical end point \24] . 

Less well studied is the bulk viscosity {() of QCD. In general, bulk viscosity vanishes 
when a system is conformally invariant such that the system is invariant under a uniform 
expansion (dilatation). For a non-interacting non-relativistic or ultrarelativistic system 
(assuming the interaction is turned off after thermal equilibrium), the system is confor- 
mally invariant and hence has zero bulk viscosity. When the interaction is turned on, 
conformal symmetry could be broken to give a finite bulk viscosity. (A notable exception 



is the infinite scattering length hmit where conformal symmetry is preserved |25l. l26l. l27l|.) 
In QCD with heavy quarks integrated out and with the hght quark masses set to zero, 
conformal symmetry is broken in the quantum level. In the perturbative region, u p t o 



some logarithmic corrections, C/s oc {1/3 



cx |28| while r]/s oc a^'^ [29, 



3Q|. 



Thus, ( is smaller than rj in the perturbative regime. When the temperature is reduced 



Tj/s reaches its minimum near T^, while (/s rises sharply near [3l|, |32, l33|. It will be 
interesting to see whether the maximum of C/s is also reached near Tc from below, which 
is the main purpose of this work. We will focus on the case with two flavors of massless 



quarks such that below Tc the dominant degrees of freedom are massless pions. 



II. LINEARIZED BOLTZMANN EQUATION AND THE GENERALIZED 
KINEMATIC THEORY 



The bulk viscosity of a system is defined by the Kubo formula 

C = - hm - / dt rfVe-* ([r;(x),T;(o)]) 



(1) 



with the trace of the energy momentum tensor. The Kubo formula involves an infinite 
number of diagrams at the leading order (LO) even in the weakly-coupled 0^ theory 
3J]. However, it is proven that the summation of LO diagrams in a weakly coupled 



theory 



34l | or in hot QED 35|, |36|] is equivalent to solving the linearized Boltzmann 



equation with temperature-dependent particle masses and scattering amplitudes. Since 
the proofs do not use properties restricted to scalar theories, the conclusion is expected 



to hold 



regime 



br mor e g eneral theories with weak couplings, including QCD in the perturbative 



28 



29 



301 1 . Here, we assume the equivalence between the Kubo formula and the 



Boltzmann equation also applies to massless pions. 

The Boltzmann equation describes the evolution of the isospin averaged pion distribu- 
tion function / = /(x, p,t) = fp{x) (a function of space, time and momentum) 



^dM^) = ^ / dTu;3p {/l/2(l + /3)(1 + fp) 
^ J 123, 



;i + /i)(i + /2)/3/p} 



(2) 



where Ep = ^/pM-mf, P = |p| and g-j^ = 3 is the degeneracy factor for three pions , 



dr 



12;3p 



1 ^ 

2Ej ' J-J- 



d^k, 



(3) 



and where T is the scattering amphtude for particles with momenta 1, 2 — > 3,p. In chiral 
perturbation theory (xPT), which is a low-energy effective field theory of QCD, the LO 
isospin averaged tttt scattering amplitude in terms of Mandelstam variables (s, t, and u) 
is 

\T\' = ^A9s' + 3{t-ur} , (4) 

where = 88.3 MeV is the pion decay constant in the chiral hmit. The pions remain 
massless below Tc and the temperature dependence of the scattering amplitude is of higher 
order and will be neglected. 

In local thermal equilibrium, the distribution function fp^\x) = ^e^(^)^f^(^)p'' — 1^ 
where (3{x) = 1/T[x) is the inverse temperature and V^{x) is the four velocity of the fluid 
at the space-time point x. A small deviation of fp from local equilibrium is parametrized 
as 



ax) = /f (a;)+5/,(a;) , 
5 ax) = -/f + (x)]xp(x). (5) 



In kinetic theory, the energy momentum tensor in a weakly interacting system is 



It is the sum of the energy momentum tensor of each particle with inter-particle interac- 
tions neglected. This is usually a good approximation when the interparticle spacing is 
much larger than the range of interaction such that the potential energy is negligible. 

The conservation of energy momentum tensor, d'^T^i, — 0, is automatically satisfied by 
the Boltzmann equation. We will decompose T^^, as 

T,, = TfJ + , (7) 

where ST^i, is the deviation from the thermal equilibrium part TjiJ . 

Tt^J = (e + V,V, - Vg,. , (8) 

where e is the energy density and V is the pressure. 

We will work at the V(x) = frame for the point x. This imphes d^y^ — after 
taking a derivative on V^j,{x)V^{x) — 1. The conservation law at local thermal equilibrium, 
d^T^^J = 0, implies 

ate+(e + P)V- V = , 
dtY + {e + V)~^VV = . (9) 



Then using the thermal dynamic relation 

e + V = T^, (10) 

one has 

(3dtY - V/3 = , 
dt/3 - /SvjV ■ V = , (11) 

where = dV/de is the speed of sound. 

The shear and bulk viscosity are defined by the small deviation away from equilibrium: 

= -2r^ (^Y^}M±1j}M _ l^^^.v . V(x)) - C%V ■ V(a:) , (12) 

where i and j are spacial indexes and Eq.l ITTi) is used to replace the time derivatives dt(3 and 
9tV by spacial derivatives V ■ V and V/3. Also, 5Toj(x) = 0, since the momentum density 
at point X is zero in the V(x) = frame. Furthermore, if there is no viscosity, the energy 
density at the same point will only be a function of T governed by thermodynamics, which 
implies STqo = 0. Viscosity could generate heat during the perturbation. However, the 
amount of heat generated should be time reversal even, because heat will be generated 
no matter whether the system is expanding or contracting. However, there is no first 
derivative term which is even under time reversal. Thus, at this order, 

SToo = . (13) 

It is easy to see why C — for ultrarelativistic and monatomic non-relativistic systems 
based on Eqs. ([6]) and f|T3|) . For ultrarelativistic systems, ~ 0; therefore, ~ by 
Eq. (j6]). For non-relativistic systems, if the particle number for each species is conserved, 
then = 26Tq = and, hence, C = 0- These are general results of the kinetic theory 
which assumes the potential energy from short-range interactions is negligible in a dilute 
system. They can be traced back to the conformal symmetry of non-interacting ultrarela- 
tivistic and non-relativistic systems. When interactions are turned on and the conformal 
symmetry is broken, Eq.([6]) has to be modified to include the effect of interaction in order 
to give the leading non- vanishing ( result. 

For pions in the chiral limit, they always satisfy the dispersion relation = even 
at finite T. This is because their goldstone boson nature prevents them from generating 



thermal masses. However, this does not imply that the system is traceless. Direct corapu- 
tation using xPT shows that trace anomaly first appears at the order of three loops 37l |. 
This is the manifestation of the gluon trace anomaly operator of QCD. In the expression 
of Eq.(l6]), = once = 0. Thus, it needs to be generalized to have non-zero T^. 
In principle, one could add two-pion, three-pion... distribution amplitudes to take into 
account the pion interaction associate with the loop diagrams. However, one can integrate 
out the medium effect and sum up the effective one-pion contributions to T^u 



T, 



TT,; 



T, 



(14) 



where T^,y is the energy momentum operator. Note that Eq. is just the leading order 
effect of the above equation which takes into account the free pion contribution to T^^ 
only. Using symmetries, T^^ has the general form: 



9-n 



(27r)3 Ep 



p^p^{i + gi{x)) + + 



(15) 



(5{xy (5{xy 

Here Lorentz symmetry is broken down to 0(3) symmetry by the temperature, and gi^^ 
are dimensionless functions of [3{x) and f-,,. In xPT, (7i_3 = (9(T'^/(47r/^)^) (37|. The 
structure {p^V^ + V^p^) is not allowed because the Tr"*" and vr^ matrix elements should 

(0)m 



be the same by charge conjugation or isospin symmetry. Thus, i^iip) 
be invariant under crossing symmetry ( p^^ —p^). In equilibrium, Tf 

e-3V 



c = Ag2 + gs 



9 



TTi{p)j should 
= e — 3V and 

(16) 



Note that energy momentum conservation is not a problem with the new terms in 
Eq. (fT5l) . In equilibrium, one just has to replace in Eq. lfTT]) by the new value to obtain 



(0) 



0. Away from equilibrium, the net effect of ( is to replace V ^ V — CV ■ V 



in Eq.([9]) which will induce second spacial derivative terms in Eq. ffTT]) . Thus, as long as 
Eq. (fT5|) gives the correct T^y, energy momentum conservation can be satisfied. 

Working to the first order in a derivative expansion, Xpi^) can be parametrized as 



Xp{x) = P{x)A{p)V-V{x)+(3{x)B{p) [pip, - -6, 



2 3 



(17) 

where A and B are functions of x and p. But we have suppressed the x dependence. 
Substituting (fTTl) into the Boltzmann equation and using Eq. ( ITTj) . one obtains one 



linearized equation for A (associated with the V ■ V structure): 

\p'-v',E'^ = SjSl [ rfri2;3p(l + ^i)(l + n2)n3(l + n,)-i 

-J ^ ^123 

X [Aip) + Aih) - Aih) - Aih)] , (18) 

where at point x, fi^\x) is written as rii = {e^^^ — l) ^ . There is also a linearized 
equation for B (associated with the (VjVj- + V jVi — trace) structure) that is related to 

;he shear viscosity 77. The computation of rj of the pion gas has been discussed in Ref. 

2o| . We will focus on solving C, in this work. 



III. VARIATIONAL CALCULATION 



Equation ([T5I 
are constants 



only determines A{p) up to a combination oi + a2-Ep, where Oi and 02 
34| . These "zero modes" (ai and a2-E'p) only appear in the analysis of bulk 
viscosity but not shear viscosity. We will discuss their effects in this section. 
The variation of Eq. (|T5|l yields 



5T, 



d^ 



^1V 5'vr 

+ u 



p 



(2vr)3Ep 



5U 



Vi^Vv (1 + ^1) + + ^2 



(19) 



f3^ f3^ 

Note that (71-3 represent loop corrections of the energy momentum tensor, thus they are 
functional of fp. To compute (, we need 



d^p 



p'{l + 9i] 



(2TifEp 

This can be simplified using the constraint 



3^2 



+ fp 



P^Sgi - 



35g2 



= 6T00 = g^ 



d^p 



(27r)3Ep 



2/11 \ , 92 + 93 
P (1+^1)+ ^2 



S92 + Sgs 



After eliminating the g2/ 0^ term in STa using the constraint, we have 



5Tii = g 

+ fp 



d3 



P 



(27r)3Ep 



4p (1 + gi) — ■ — 
g2 + g3 



p^^gi 



45-2 + 5-3 3 (5(2(55(3 - 53^52) 



+ 



g2 + 53 {.g2 + 53) 
2^4 



^9.d I ^,p5fp (1 + Oi 



(47rA)^ 



(20) 



(21) 



(22) 



where d = {4g2 + 93) / [92 + 93) and the pion remains massless in the chiral hmit even at 
finite T, so we have used = E^. The above expression for ST^ imphes 

C = ^9.(3d J ^EpTip (1 + rip) Aip) . (23) 



Then using Eq. f|T8l) and the symmetry property of the scattering amphtude, 

i=l,2,3,p 



X (1 + m) (1 + n2)n^np [A{p) + A{ks) - A{k2) - A{h)f . (24) 

Note that equating Eqs. (123!) and fl24|) is equivalent to taking a projection of Eq. (ITSjl . 
t can be shown that any ansatz satisfying Eqs. ( 123|) and ( 12^ gives a lower bound on ( 



381]. Thus, one can solve ( variationally, i.e. finding an ansatz A{p) that gives the biggest 

C- 

It is known that if one uses the ansatz A{p) = oi + a2-E'p, then it will not contribute 
to the 2^2 scattering on the right-hand side of Eq. (fT8|) (the 02 terms cancel by 
energy conservation). In fact, this ansatz will not contribute to all the particle number 
conserving processes but can contribute to particle number changing processes, such as 
2 4 scattering, which we have not shown. As we know from Eqs. ( fT8|) and ( l23l) . ( 
is proportional to the size of A{p) which is inversely proportional to rate of scattering. 
Thus, if the 2-^2 scattering has a bigger rate than the 2^4 scattering, then this ansatz 

gives a bigger ( by bypassing the faster 2-^2 scattering. In 0^ theory, it was found that 
( is indeed set by the 2^4 scattering [34]. However, in perturbative QCD (PQCD), the 



soft particle number changing bremsstrahlung is faster than the 2—^2 scattering 281]. 
Thus, ( is governed by 2 ^ 2 scattering. 

In the case with massless pions, however, 2 — 2 scattering is still the dominant process. 
While using the ansatz A{p) = ai + 02 -Ep, the 6Toq = constraint in Eq. fl2T]) demands 
O1/02 = because Up oc 1/p as p — > 0. Since A{p) parametrizes a small deviation of fp 
away from thermal equilibrium, ai/a2 = gives ai = instead of 02 ^ C)0 and ai finite. 
Thus, to maximize (, we uses the ansatz A{p) = a2Ep + asi^^ + ... without the ai term. 
The point is, 2 — > 2 scattering cannot be bypassed and it will be the dominant process in 
our calculation. 

To compute C, it is easier to eliminate the (1 + 91) term in Eq.( l20i) using Eq. (l2Ti) : 



STii 



+ fp 



A6g2 + 6gs 



5L 



4^2 + 93 



(25) 



Note that g2 and terms at 0{T^/ (4-7? /^r)^) arise from three-loop diagrams and from two- 
loop diagrams with insertions of higher order counterterms and each loop integral has one 
power of fp in the integrand. Thus, we will make an approximation here to assume the 
{^5g2 + 5g^) term is proportional to the 5fp term with a proportional constant (/ — 1), 
where I means the power of fp (or the number of loops) in Ta. Since I is between 2 and 3, 
we take the mean value / = 2.5 and associate the uncertainty of / to the error estimation 
ofC 

The trace anomaly for massless pions appears from three-loop diagrams and from two- 
loop diagrams with insertions of higher order counterterms [all are 0{T^ / f^)] 37|. Thus, 

For two-loop diagrams, the associated / factor is 2 while for three-loop diagrams, the 
associated / factor is 3. Here, without distinguishing the contribution from each diagram, 
we take the mean value / = 2.5 and associate the uncertainty of / to the error estimation 
of Thus, 

C = / ^An, (1 + Up) A{p) . (26) 



3/5 J (27r)3Ep 



Note that A{p) (x g^^ — vf) from Eq. (|T8l) . Thus, for massless pions. 



C = hl 



3^) ( 3 - v: 



J n 
2^5 



(27) 



where is given by dimensional analysis and h is a dimensionless constant. To find the 
numerical solution for h, we neglect the higher-order gi^^ terms in Eq.( !2T]) and use the 
ansatz A{p) = Yl^=i ^nV^- We find 

/i ~ 65 . (28) 



Using the xPT result of Ref. [37j for e and P, we obtain 



c--K^)('4-^)('4-^)S 

where Ap ~ 275 MeV. As expected, the bulk viscosity vanishes 
coupling between pions vanishes. 



as 



(29) 

oo or when the 



The leading order contribution for pion entropy density s is just the result for a free 
pion gas: 



45 



(30) 



The trace anomaly for massless pions appears from three-loop diagrams and from two- 
loop diagrams with insertions of higher order counterterms [all are 0{T^/ f^)] 37|. Thus, 

For two-loop diagrams, the associated / factor is 2 while for three-loop diagrams, the 
associated / factor is 3. Here, without distinguishing the contribution from each diagram, 
we take the mean value / = 2.5 and associate the uncertainty of / to the error estimation 
ofC 



ST ~ 
"-^ /li/ — 



gJ 



d^p 



(27r)3Ep 



P^Py (1 + ^l) + 



Using Eqs. ([TB]) and ([T^, 

5Tiiix) 



STJx) 



gn 



= -gn 

+ u 

d^p 

(27r)3Ep 



d3 



p 



(27r)3Ep 



5L 



/32 

45^2 + gs 

Pixy 



(3^ 



{5U[Ap\l+g,){Ag2 + g, 



(45-2 + 5-3) p^^gi + 



3 {g25gi - gi5g2 
P^ 



(31) 
(32) 

(33) 
(34) 



where the factor / appears because in the computation of STu from Eq. f|T5|) . not only 
fp but also (71-3 depend on 5fp. In xPT, the trace anomaly for massless pions appears 
from three-loop diagrams and from two-loop diagrams with insertions of higher order 
counterterms [all are 0{T^ / f^)] [37|. For two-loop diagrams, the associated / factor is 2 
while for three-loop diagrams, the associated I factor is 3. Here, without distinguishing 
the contribution from each diagram, we take the mean value / = 2.5 and associate the 
uncertainty of / to the error estimation of C,. 

The above-mentioned constraint, 5Tqq = 0, yields. 

The dimensionless combination C/'^ is shown in Fig. 1. The solid line below Tc is the 
leading order massless pion gas result (we have used / = 2.5, explained below Eq. (12CT]) . 
and the lattice result, ^ 200 MeV, for 2+1 flavors of improved staggered fermion as an 
estimation 39j). The error on this curve is estimated to be 30%-40% from / and higher- 



order corrections. But the monotonic increasing behavior should be robust. The solid 
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FIG. 1: (Color online) C,/ s shown as a function of T jT^. The solid line below Tf. is the massless 
pion gas result (Tf. ~ 200 MeV and I = 2.5, explained below Eq. (j20p . are used). The error on 
this curve is estimated to be 30%-40%. The points are the lattice results for gluon plasma 2]. 
The solid and dashed lines above Tc give the central values and the error band from the QGP 



sum rule result of Ref. 



22]. 



points are the lattice results for gluon plasma [33;]. The QGP curves above Tc (the solid 
line gives the central values and the dashed lines give the estimated errors) are based on 
an exact sum rule, a lattice result for the equation of state, and a spectral function ansatz 
with massive quarks 3^. Since the light quark mass dependence in the QGP curve is 
expected to be small. Fig. 1 shows that, in the chiral limit, QCD C,/ s reaches its maximum 
while T]/ s reaches its minimum around Tc as mentioned above. The same C,/s behavior is 



40| 



also seen in molecular-dynamics simulations of Lennard- Jones model fluids 

A recent massive pion gas calculation shows that C, has two peaks 43|, one is near 
10 MeV and the other is near Tc. They are corresponding to breaking of the conformal 
symmetry by the pion mass and the anomaly, respectively. The behavior near the higher 
temperature peak is similar to what we have found here for the massless pion case. It is 
also similar to the C, behavior of near Tc with Hagedorn states included. The behavior 
near the lower temperature peak is similar to earlier results of 
calculation of 43[ also conforms our qualitative behavior of C,. 



45 



46j. The massless pion 



In the large Nc (the number of colors) limit, 



and 



C 1 

- oc p^2]\f2 massless pion gas, (35) 
^ f 



^ cx ^ cx for PQCD, (36) 



where we have used the scaling f^, oc a/ZVc, gn oc Nj, oc l/N^ and Nf is the number of 
light quark flavors. Also, for massless pions. 



^-i**n^ju-7ju-"-j • '''' 

This is similar to C/v ^ 15(1/3 — f^)^, which is obtained for a photon gas coupled to hot 
matter 411] and is also parametrically correct for PQCD [28]. This is because in those 
cases, 2 — i> 2 scattering is the dominant process in both ( and t] computations. It is not 
the case, however, in (f)^ theory in which (1/3 — ff) C/v has large T dependence because 
( is dominated by 2 4 scattering while t] is dominated by 2 — > 2 scattering. The scaling 
is also different from (/rj oc (1/3 — f^) for strongly coupled J\f = 2* gauge theory using 
AdS/CFT [421. 



IV. CONCLUSIONS 



We have computed the bulk viscosity for a gas of massless pions using the Boltzmann 
equation with the kinetic theory generalized to incorporate the trace anomaly. The re- 
sulting (Js, together with the corresponding results of gluon plasma [3^ and quark gluon 



plasma 3l| indicates (/s reaches its maximum near Tc while rj/s reaches its minimum 



near T^. If the (/s behavior is unchanged for massive pions, then the hadronization of the 



fire ball in heavy ion collisions would imply large entropy production [31|, |33j and slow 
equilibration. It would be interesting to explore the implications of the possible large bulk 
viscosity near a phase transition in cosmology if the phase transition above the TeV scale 
is based on some strongly interacting mechanism. 
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